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Abstract 

Reaction-diffusion equations based on a poly- 
merization model are solved to simulate the 
spreading of hypothetic left and right handed 
life forms on the Earth's surface. The equa- 
tions exhibit front-like behavior as is famil- 
iar from the theory of the spreading of epi- 
demics. It is shown that the relevant time 
scale for achieving global homochirality is not, 
however, the time scale of front propagation, 
but the much longer global diffusion time. The 
process can be sped up by turbulence and 
large scale flows. It is speculated that, if the 
deep layers of the early ocean were sufficiently 
quiescent, there may have been the possibil- 
ity of competing early life forms with oppo- 
site handedness. Key Words: Homochirality 
Origin of life Exobiology. 

Introduction 

There has been a remarkable development in assess- 
ing the probability of extraterrestrial life in the Uni- 
verse. Although higher life forms are now believed 
to be exceedingly rare, primitive life forms may be 
more wide spread than what has been believed some 
20 years ago (Ward and Brownlee, 2000). Two major 
discoveries have contributed to this view: (i) micro- 
fossil and carbon isotope evidence that life has ex- 
isted on Earth at least as early as 3.8 billion years 
ago (Schopf, 1993; Mojzsis et ai, 1996), and (ii) the 
discovery of extremophilic life forms on Earth making 
it plausible that life can exist (or has existed) on other 



celestial bodies in our solar system and beyond. How- 
ever, Cleaves and Chalmers (2004) present a more 
cautions view quoting evidence that some of the old- 
est viable organisms were actually only mesophiles, 
and not at least hyperthcrmophiles. 

Life may either have originated elsewhere and 
then delivered to Earth [but again, see Cleaves and 
Chalmers (2004) for a more cautions view] , or it may 
have been created locally from prebiotic chemistry, 
e.g. in hydrothermal systems at the ocean floor (Rus- 
sell and Hall, 1997; Martin and Russell, 2003). In 
either case there is the possibility of multiple sites 
on Earth, or alternatively in space, where different 
life forms may have coexisted for some time - even 
though on Earth only one life form has prevailed in 
the end. Indeed, the DNA of the most primitive life 
forms on Earth suggests that it is based on a sin- 
gle common ancestor (e.g., Collins et ai, 2003). One 
important characteristics of life, as we know it, is 
the handedness (or chirality) of amino acids and sug- 
ars. It is well known that the sugars generated in 
living organisms are dextrorotatory. Examples in- 
clude the glucose produced in plants via photosyn- 
thesis, but also the ribose sugars in the backbones 
of the nucleotides of the DNA and RNA. Conversely, 
the amino acids in proteins of living organisms are all 
left-handed. 

A commonly discussed pathway for the origin of life 
is a pre-RNA world (Bada, 1995; Nelson et ai, 2000), 
that could have been - in its earliest form - achi- 
ral. An early peptide world composed of glycines has 
been discussed (Milner- White and Russell, 2004) and 
also peptide nucleic acids (PNA) based on a glycine 
backbone are achiral (e.g., Pooga et ai, 2001). PNA 



molecules are frequently discussed for transcribing 
the genetic code of early life forms (Nielsen et al, 
1993) and even artificial life forms (Rasmussen et al., 
2003). However, PNA based on other amino acids 
is chiral (Tedeschi et al, 2002), and it might there- 
fore be plausible that the onset of chirality could have 
emerged during the evolution of the PNA world and 
hence after the onset of life. 

There is a priori no reason why life on Earth 
could not have been is based on nucleotides with left- 
handed sugars that catalyze the production of right- 
handed amino acids, i.e. just opposite to what it actu- 
ally is. It is now widely accepted that at some point 
around the time when the first life has emerged (or 
some time thereafter) there must have been a random 
selection favoring one particular chirality. Whether 
it was right or left handed must have been a matter 
of chance and can be explained by self-catalytic poly- 
merization of nucleotides (Frank, 1953). An impor- 
tant aspect of this process is the cross-inhibition by 
nucleotides with the opposite chirality which would 
spoil the polymer and prevent further polymeriza- 
tion on the corresponding end of the polymer (Joyce 
et al., 1984). Although polymerization of alternating 
left and right handed nucleotides is usually impos- 
sible (Kozlov et al, 1998), this can be different for 
amino acids in more complicated molecules, for ex- 
ample certain siderophores ( "iron bearers" ) can con- 
tain mixed left and right handed amino acids (Mar- 
tinez et al., 2000). Another example are the so-called 
"nests" , i.e. are pairs of left and right handed amino 
acids within relatively short sequences (Watson and 
Milner-White, 2002; Pal et al, 2002). 

The first detailed polymerization model taking 
such cross-inhibition into account was proposed by 
Sandars (2003). This model was studied further by 
Brandenburg et al. (2004a, hereafter referred to as 
BAHN) and Wattis and Coveney (2004). Since the 
influential paper of Frank (1953) the various models 
that have been devised to explain the origin of the 
same chirality (or homochirality) of nucleotides and 
amino acids were all based on a pair of two ordinary 
differential equations describing the time evolution of 
some collective properties of the two types of poly- 
mers with left and right handed building blocks (e.g., 
Saito and Hyuga, 2004a; hereafter referred to as SH). 



As was shown in the paper of BAHN, a modified ver- 
sion of the model of SH can actually be derived as a 
reduction of the full polymerization model of Sandars 
(2003). 

In the present paper we extend the model of BAHN 
by including the fact that opposite chiralities are 
likely to originate at different places. At those places 
where life forms of opposite chirality meet, one of 
the two will eventually dominate and extinguish the 
other. If the time for this to happen was sufficiently 
long, life forms of opposite handedness may have co- 
existed for some time at different places on the early 
Earth. 

Mathematically speaking, our model falls in the 
class of diffusion-reaction (or advection-diffusion- 
reaction) equations that arc being discussed in many 
different fields; see the book by Murray (1993) for bi- 
ological applications. In this paper the reactions cor- 
respond to the auto-catalytic polymerization which 
corresponds to a local instability with a given linear 
growth rate. The importance of spatial extent has al- 
ready been emphasized by Saito and Hyuga (2004b) 
who generalized the SH model by using a Monte- 
Carlo method as it is used in percolation studies. 

Wc begin by outlining the basic alterations to mod- 
els without spatial extent and discuss relevant time 
scales for achieving homochirality in a spatially ex- 
tended system. We then briefly review the polymer- 
ization model of Sandars (2003) and the reduction 
of BAHN, discuss the effects of spatial extent and 
present results of simulations. 

Relevant time scales 

Since the seminal paper by Frank (1953) we know 

that, in the case without spatial extent, an infinites- 
imally small initial enantiomeric excess rj (defined in 
the range — 1 < < 1; see below) can grow exponen- 
tially like 

?7 = ?7oe^* (for |77-77o| « !)• (1) 

Here ?7o is the initial value, which must be close to 
the racemic value of zero (see BAHN for the calcula- 
tion of A). The growth rate depends on the reaction 
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coefficients and the source term replenishing the sub- 
strate, from which new monomers can grow. In the 
case with spatial extent, however, there is the pos- 
sibility that life forms with opposite handedness can 
emerge independently at different sites. This can sub- 
stantially delay the evolution toward homochirality, 
because now the relevant time scale is governed by 
the time it takes to extinguish already existing life 
forms of opposite handedness. 

The situation is somewhat analogous to the propa- 
gation of epidemics whose expansion over the Earth's 
surface is known to be related to the speed at which 
fronts propagate in reaction-diffusion systems. A pro- 
totype of such a system is the Fisher equation (Fisher, 
1937; KaUappan, 1984; Murray, 1993) 

where c is the concentration, A is the growth rate in 
the homogeneous limit (c independent of a;), cq is a 
saturation amplitude, and k is a diffusion coefficient. 
Equation |(2J) admits propagating wave solutions of 
the form c = c{x — Wfront^) where 

t;front ~ 2(«A)i/2 (3) 

is the speed of the front [Murray (1993); see also 
Mendez et al. (2003) for the case of a heterogeneous 
medium] . A typical application of reaction-diffusion 
equations is the spreading of the Black Death in Eu- 
rope in the fourteenth century. In this connection, 
Murray (1993) quotes an effective diffusion coefficient 
of about 3 X IC^km^/yr (corresponding to IQ-^m^/s) 
and a reaction rate (or in this case a 'mortality rate') 
of 15 yr~^ (corresponding to A ~ 5 x 10^^ s~^), which 
gives Wfront ~ 5cm/s or about 1300km/yr. 

Applying the same front speed (1300km/yr) to the 
homochiralization over the whole Earth surface (typ- 
ical scale L ^ 10 000 km), one finds a typical time 
scale Tfront — i/f front, which would be on the order 
of 10 years. However, it will turn out that the time 
scale involving front propagation is too optimistic an 
estimate for the problem at hand, because once right 
and left handed life forms meet, the reaction front 
will almost come to a halt. As it turns out, the only 
reason why these fronts can still propagate is because 



of their curvature. This leads eventually to the com- 
plete dominance of only one handedness. The time 
scale for this to happen turns out to be the global 
diffusion time scale 

Tdiff = L^/n- (4) 

Putting in the same numbers as above one finds ~ 
3000 yr. This is still rather short, but this is because 
the adopted value for k is more like an effective (eddy) 
diffusivity and hence rather large compared with the 
usual molecular values quoted by Cotterill (2002) for 
air (k = lO^^m^/s) or water {k = 10~^m^/s). These 
values would lead to diffusion times in excess of the 
age of the universe. 

Larger values of n are probably more realistic in 
that they take into account the effects of additional 
macroscopic mixing processes that are required to es- 
tablish global homochirality. Examples include ocean 
currents, surface waves, as well as winds and turbu- 
lence in the atmospheric boundary layer (e.g. Garrett, 
2003) . Below we shall refer to the corresponding time 
scale as the turbulent time scale 

TtvLvh = L/Ur^s, (5) 

where Urms = {u^)^^^ is the root mean square veloc- 
ity. 

A commonly used procedure to estimate turbu- 
lent diffusion times is in terms of enhanced turbu- 
lent (eddy) diffusivities. We return to this in more 
detail, but for now we just note that in the ocean 
values between ~ lO^^'m^/s for the vertical eddy 
diffusivity and up to ~ lO'^m^/s for the horizontal 
eddy diffusivity are typical (see Chap. 8.5 of Stewart, 
2003). This would again suggest turbulent diffusion 
times of around 3000 years, which would then be the 
time required for achieving global homochirality in 
the ocean. 

The estimates outlined above can be put on a more 
quantitative basis. Before we can study the analogous 
equation to the Fisher equation, we briefly review the 
basis of the polymerization model of Sandars (2003) 
and its reductions as proposed by BAHN. Since we 
are interested in the spatially extended case, we add 
the effects of diffusion and advection to these equa- 
tions. We therefore consider the case where the con- 
centrations of the different polymers [Ln] and [i?n]. 



3 



obey advection-difFusion-rcaction equations, i.e. ev- 
erywhere in the equations in the paper by BAHN we 
replace 



(6) 



Hero, u and k arc the velocity vector and the diffusion 
coefficient, respectively. 

The polymerization model 

In the polymerization model of Sandars (2003) it is 
assumed that a polymer L„ with n left-handed build- 
ing blocks can grow on either end by a left or right 
handed monomer, Li or Ri, respectively, with a re- 
action coefficient 2ks or 2fc/, respectively. The latter 
rc;action is referred to as enantiomeric cross-inhibition 
leading to a polymer whose one end is now 'spoiled', 
allowing further growth only on the remaining other 
end at rates ks or ki for a reaction with Li or 
respectively. (The left handed polymers L„ are not 
to be confused with the size of the domain L.) The 
full set of reactions included in the model of Sandars 
(2003) is thus 



Ln + Li 


2fes 


Ln+1, 


(7) 


Ln + Rl 


2ki 


LnRl, 


(8) 


Ll + LnRl 


ks ^ 


Ln+lRl, 


(9) 


-Rl + LnRl 


ki ^ 


RlLnRl, 


(10) 



where ks and kj are the reaction coefficients for 
monomers to a polymer of the same or of the op- 
posite handedness, respectively. For all four equa- 
tions we have the complementary reactions obtained 
by exchanging L ^ R. 

In addition, new monomers of either chirality are 
being regenerated from a substrate 5 at a rate that 
may depend on the already existing excess of one chi- 
rality (i.e. the enantiomeric excess), so 



Ru 



S"'^ Li, 



(11) 



where kc is proportional to the growth rate of 
monomers and Cr and Cl determine the enzymatic 



enhancement of right and left handed monomers. 
Sandars (2003) assumed that Cr = Rn and Cl = 
Ln, but different alternative proposals have been 
made, for example = ^ i?„ and Cl — Ln 
(Wattis and Coveney 2004) or Cr = Yli^^n and 
Cl = Y^nLn (BAHN). The overall behavior of 
the model is not however affected by these different 
choices. 

The governing equations, ignoring spatial depen- 
dencies, are (for n > 3) 

'^^=2ks[L{\{[Rn-i]-[Ln])-2ki[Ln][Ri], (12) 

^ = 2ks[Rl]{[Rn-l] - [Rn]) - 2ki[Rn][Li], 

(13) 

where the factor 2 in front of ks and ki reflects the 
fact that the monomers can attach to either side 
of the polymer. For n = 2, however, to identical 
monomers interact, and the number of correspond- 
ing pairs is not m^, as for non- identical pairs (n > 2), 



but only ^m{in — 1) 



1 2 



canceling the factor 2 



in front of the [Li]^ term (see Sandars 2003, and also 
BAHN). Thus, we have 



d[L2] 
dt 

d[i?2] 
dt 



ks[Li]^ - 2fcs[Li][L2] - 2ki[L2][Ri], (14) 



= ks[Ri]'' - 2ks[Ri][R2] - 2ki[R2][Li]. (15) 



The evolution equations for [Li] and [Ri] are given 

by 



d[Li] 
dt 

where 



Ql-Xl[Li], 



d[Ri 
dt 



Qr-Xr[Ri], (16) 



N-l 



N 



XL = 2ksJ2i^n] + 2klY,[Rn] 
n=l n=l 
N-l N 
+ksJ2[LnR] + kjJ2[RnL], 



(17) 



n=2 



n=2 



and 



N 



XR = 2ksJ2[Rn]+2ki^[Ln] 



n=l 



4 



N-1 



N 



+ ks J2[RnL] + kiJ2[LnR], 



(18) 



n=2 



Tl = 2 



are the decay rates that quantify the losses associated 
with the reactions f7jl- Hl()(l . respectively. Here N = 
iV - 1 in the model of Sandars (2003) and AT = TV in 
the model of BAHN. 

New monomers are generated from a substrate 
with concentration [S] via the source terms Ql and 



}L^kc[S]ipCL+qCR), 
^R^kc[S]ipCR + qCL), 



(19) 
(20) 



where kc is the corresponding reaction coefficient, 
p — ^(1 + /) and q = i(l — /), and / is the fidehty 
of the enzymatic reactions. The concentration of the 
substrate, in turn, is maintained by a source Q, via 



d[S] 
dt 



Q-IQl + Qr), 



(21) 



where Ql + Qi? = kc[S]{CL + Cr); see Eq. (EOJ. In 
general, we expect Cl and Cr to be functions of L„ 
and Rn, respectively. Sandars (2003) assumed Cl = 
[Ln] and Cr = [Rn], i-e. the catalytic effect depends 
on the concentrations of the longest possible chains 
of left and right handed polymers, but other choices 
are possible (see BAHN) . The degree of handedness 
(or enantiomeric excess) is generally characterized by 
the ratio of the difference to the sum of the number 
of left and right handed building blocks, i.e. 



V = 



Er — El 
Er + El 



(22) 



where Er = ^7i[i?„] and El = 'Y^n[Ln] measure 
the total number of left and right handed building 
blocks. 

The important point is that the polymerization 
equations show bifurcation behavior where / is the 
control parameter and rj is the control parameter. 
There is a critical value, /ditj above which the 
racemic solution (77 — 0) is linearly unstable (with 
a growth rate dlnrj/At = A) and a new solution with 
enantiomeric excess is attained. Only for / = 1 is 
the enantiomeric excess complete, i.e. 77 = 1 or —1. 



The value of /crit depends on the ratio kj/ks- For 
example in the model of BAHN, using kj/ks = 1, we 
have /crit ~ 0.38, and /crit increases for decreasing 
kj/ks- For kj/ks = (no enantiomeric excess) no 
bifurcation is possible, i.e. /crit = 1- 

The polymerization problem without spatial extent 
is governed by the following characteristic time and 
length scales 



{Qks) 



-1/2 



{Q/ksV\ (23) 



where the subscript 'micro' is chosen to emphasize 
that with spatial extent the problem will be governed 
by other 'macroscopic' time scales such as Tdig and 
Tturb (see previous section). 

Since neither Q nor ks are well known, it may be 
advantageous to assume Tmicro and ^micro instead: 



micro micro' 



^ ^micro mici 



(24) 



A plausible time scale for polymerization is hours 
(Tmicro ~ 10^ s. Say), but it could be much longer, and 
typical molecule concentrations are 10^° cm""^, but it 
could also be much less. The nominal values corre- 
sponds to ks ^ 10~^^ cm^ s~^, Q ~ 10^^ cm~^ s~^. 

The growth rate of instability of the racemic so- 
lution, A, scales with t^^^^, and it is therefore con- 
venient to introduce the symbol Aq = T^lcm- ^^^^ 
example, for kj/ks = / = 1, we have A = 0.57Ao; see 
BAHN. 

Allowing for spatial extent 

Before discussing the reduced model and its extension 
to allow for spatial extent, let us briefly illustrate the 
principal effects of spatial diffusion on the full poly- 
merization model. If the fidelity is close to unity, an 
initially small perturbation to the racemic solution 
([i?„] = [L„] = 2-"[Li] for ki = ks and n > 2; see 
BAHN) grows until a new homochiral state is reached 
where there are polymers of only one handedness with 
equally many of each length (i.e. [L„] = 2[Li] for 
n > 2). The same can happen also locally if we al- 
low for spatial extent provided the domain is large 
{L ^ Cmicro)- In Fig. ^we show the result of a one- 
dimensional numerical calculation in a periodic do- 
main — L/2 < X < L/2 with £/^micro = 10 and finite 
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plots of [Rn 



and 

[Ln] for t/rdiff = 0.8 as a function of x 
and n, and the corresponding dependencies of 
I]^=i'^[^"] (solid line), compared 

with \NY.'^^^[Rn] and iiV^^^^^ii,,] (dashed line), 
and \n'^[Ri] and |iV2[Li] (dotted line), aU in 
units of {Q/ksY^^ ■ The normalized diffusivity is 
k/{L^Xo) = 10-2 and N = 20. 



diffusivity k/{L'^Xq) = 10 2, i.e. Tdiff = lOOrmicro, 
but no flow {u — 0). 

Note that the domain divides into regions with 
opposite handedness. If one waits a few e-folding 
times (A~^), each of the homochiral sub-domains pos- 
sesses polymers of appreciable length (Fig. At 
later times the polymers will have grown to maxi- 
mum length, with a sharp interface separating do- 
mains with right and left handed polymers. The evo- 
lution of these domain walls will be discussed in de- 
tail in the next section using the reduced model of 



BAHN. 



The reduced model 

It turns out that, without changing the basic proper- 
ties of the model, a minimal version is still mean- 
ingful for N = 2, and that the explicit evolution 
equations for the semi-spoiled polymers, [L2R1] and 
[R2L1], can be ignored. However, the associated re- 
action (and likewise for L ^ R) still leads to loss 
terms in all equations (proportional to kj). In fact, 
the very presence of the kj terms is essential for the 
bifurcation toward homochirality. 

Thus, we only solve Eqs. 1)141) and (|15|l together 
with Eqs. JTSJ-lIHIl- Following Sandars (2003), we 
also assume that Cl = [L2] and Cu ~ [R2] (in- 
stead of Cl — El and C_r — Ejf , which would yield 
more complicated expressions). A further simplifi- 
cation can be made by regarding [L2] as a rapidly 
adjusting variable that is enslaved to [Li] (and sim- 
ilarly for [i?2])- This technique is also known as the 
adiabatic elimination of rapidly adjusting variables 
(e.g., Haken, 1983). Equations (|12l) and H13|) become 

= ks[Li]^ - 2[L2](fcs[ii] (25) 
= fcs[i?i]' - 2[i?2] (ks[Ri] + , (26) 



which are solved for [L2] and [R2], respectively. These 
quantities couple back to the equations for [Li] and 
via Ql and Qli.. Finally, we also treat the sub- 
strate [S] as a rapidly adjusting variable, i.e. we have 
kc[S] = Q/{[L2] + [i?2])- We emphasize that the 
adiabatic elimination does not affect the accuracy of 
steady solutions. It is convenient to introduce new 
dimensionless variables, 

X=[i?l](2fc5/g)^/^ Y=[L,]{2ks/Q)'^^. (27) 

and a constant Aq — {Qks/2y^^ , which has dimen- 
sions of (time)-^. The introduction of the l/\/2 fac- 
tor is useful for getting rid of factors of 2 in the re- 
duced model equations below. 

In order to compare first with the SH model we 
restrict ourselves to the special case kj/ks = 1, which 
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leads to the revised model equations 

= E^l±l^-rX. (28) 

where dots denote derivatives with respect to time 
and r = X + Y, = + Y^, p= {1 + f)/2 and 
(7 = (1 — /)/2 have been introduced for brevity. The 
reduced equations resemble those of the SH model 
in that both have a quadratic term proportional to 
X'^ (or Y'^), which is quenched either by a 1 — r fac- 
tor (in the SH model) or by a 1/{X'^ + F^) factor in 
our model. Furthermore, both models have a back- 
reaction term proportional to —X (or —Y), but the 
coefficient in front of this term (denoted by A in the 
SH model) is not constant but equal to {X + Y). 

The form of the equations 1)28(1 and H29|l suggests 
a transformation to a new set of variables 



S = X + Y, A = X -Y, 



(30) 



in terms of which the enantiomeric excess is in the 
reduced model simply 



■n = A/S. 



(31) 



Adding and subtracting Eqs. ((28|l and ((29|l from each 
other leads to a new pair of equations in terms of S(t) 
and A{t): 



AqU = 2/ 



SA 



52 + ^^2 



-SA. 



(32) 
(33) 



Fixed points are: S = 1, A = 0, ±\/2f — 1, cor- 
responding to the points ((1 ± \/2/ — l)/2, (1 =F 
V2/- l)/2), (i, i) in the (X,y)-plane (the S = -1 
fixed points are not physical, because < X,Y < 1). 

We see that the enantiomeric excess is complete 
only for / = 1 and that for this case physical sta- 
ble points with an enantiomeric excess only appear 
with / > 1/2 (for kj/ks = 1). The growth rate of 
the enantiomeric excess, 77, from the racemic point, 
{S, A) — (I, 0), can be found by linearizing the equa- 
tions and it turns out to be ?) = 2/ — 1. For simplicity, 
we concentrate here on the limiting case / = 1, but 



also show in a later section, how the results depend 
on changing / (and hence the growth rate). 

The first equation, Eq. (|33|) . is easily solvable in 
terms of hyperbolic functions, 

_ r cothAo(i - to) for 5(0) > 1, 
'^^^^ ~ \ tanh Ao(i - to) for 5(0) < 1, ^'^ > 

where to is an integration constant that can also be 
expressed in terms of the initial value 5(0). From 
Eq. H34|) we note the long time behavior: 5^1, 
regardless of the initial value. Since < X,Y < 1, 
we have < S < 2. The corresponding solution for 
A, assuming / = 1, is 



a — \ a? + cosh^ Ao(t — to) 
A{t) = ^ . , ^ (35) 



sinh Ao(t — to) 



for 5(0) > 1, and 



a — \/ + sinh^ Ao(t — to) 

At) = r^^ ^ (36) 

coshAo(t-to) ^ ' 

for 5(0) < 1. Here a is another integration constant 
that depends on the initial condition: 

- [^(0)2-5(0)2] (37) 



2^(0) 



for 5(0) > 1, and 



« = +£^!ii^ [^(0)2 _ 5(0)2] (38) 



2AiO) 



for 5(0) < 1. 



The reduced model with spatial extent 

As in the model without spatial extent, the polymer- 
ization dynamics seems to be controlled by the dy- 
namics of the shortest polymer. In other words, the 
polymerization dynamics of the longer polymers is 
enslaved by that of the shorter ones. 

This is the basis of the reduced model of BAHN 
and hence we consider a reduced model with spatial 
extent. The new set of reduced model equations reads 
as 

dt+u-V -kV^)X = XoFx{X,Y), 
dt + u-V - k\7^)Y ^ XoFy{X,Y), ^ > 
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where Fx and Fy are the right hand sides of the 
reaction equations (|28() and H29I) . 

A crucial property of this model is its bifurcation 
behavior in time, as we have seen when considering 
the models with no spatial extent. For u = k = 
0, we now have the spatially homogeneous solution 
X{x,t) = Y{x,t) — 1/2, i.e. there are equally many 
right and left handed molecules. Furthermore for / > 
1/2, there are the two additional solutions. 



where 



X{x, i) = 1 - Y{x, t)^~l± y/2f - 1. (40) 

The question we wish to answer is how an initially 
nearly homogeneous racemic mixture attains a final 
state, whether or not it is composed of regions of 
different handedness and what are the relevant time 
scales. 

The one-dimensional case 

Using the previously introduced variables S, A and 
specializing to the one-dimensional case, the equa- 
tions without advection (|39|l take the form 



S = Xo{l - S'^) + kS" 



A = XqSA 



(41) 

1 ) + kA", (42) 



where primes indicate a derivative with respect to x. 
Clearly S{x,t) = 1 is a solution of the first equa- 
tion. Numerically we have seen that S tends to 1 
very quickly from a 5 ^ 1 initial state. Unless k is 
very large, one can view the problem in two parts: 
first S relaxes to 5 = 1 and then A begins to evolve 
spatially. With 5=1, Eq. (gJl takes the form 



A = XoFA{x,t) + kA", 



where we have introduced 



Fa = A 



i^A" 
l + A^' 



(43) 



(44) 



The time independent equation A — can be inte- 
grated once and written in a suggestive form as 



V{A) = \A^ - ln(l + A^) + Vq (46) 





2Ao 



(A'f - V{A) = 0, 



(45) 



is the potential and Vb is an integration constant. 
The form of the potential is shown in the left hand 
panel of Fig. [3 The minimum is degenerate, the 
two minima are located at ^ — ±1 reflecting the 
Z2 symmetry of the theory. The minima correspond 
to the two stable solutions the system tends to. As 
usual, there is a domain wall solution interpolating 
between the two minima. The solution is shown in 
the right hand panel of Fig. 13 The shape is very well 
approximated by A(x) = tanh(— ax-\/Ao/K), where 
a ~ 0.58. So in one dimension the initially racemic 



< 
> 



Figure 2: Potential V{A) (for Vq = 0) and the solu- 
tion interpolating between the left and right handed 
regions. 

mixture tends to form domain walls separating left 
and right handed regions. One can therefore see the 
process as spontaneous symmetry breaking into dif- 
ferent domains. Numerical work confirms the analyt- 
ical considerations: the region breaks into different 
domains that once relaxed to 77 = ±1, are separated 
by time independent domain walls. In one dimension, 
it is hence difficult to achieve a homogeneous mix- 
ture unless the initial perturbation is only directed 
towards one of the minima. 

In Fig. 131 we show a space-time diagram of the solu- 
tion of the one-dimensional problem without advec- 
tion and an initial perturbation corresponding to a 
weak right-handed excess at x/L — 0.1 and a some- 
what stronger left-handed excess at x/L = —0.1. 
Note the propagation of fronts with constant speed 
Vbonti as given by Eq. Q, if the exterior is racemic 
(i.e. X = Y = 1/2) and a nonpropagating front when 
the chirality is opposite on the two sides of the front. 
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-0.4 -0.2 0.0 0.2 0.4 




Figure 3: Profile of X{x, t,) and space-time diagram 
of X(x^t) for the one-dimensional problem without 
advection and an initial perturbation corresponding 
to a weak (amplitude 0.01) right-handed excess at 
x/L = 0.1 (marked in white or yellow) and a some- 
what stronger (amplitude 0.3) left-handed excess at 
x/L = —0.1 (marked in dark or blue). Note the prop- 
agation of fronts with constant speed if the exterior is 
racemic (i.e. X = Y = 1/2, shown in medium shades 
or red) and a nonpropagating front when the chirality 
is opposite on the two sides of the front. The normal- 



ized diffusivity is k/(L^Ao) 
in Fig. [2 



10- 



i.e. the same as 



The higher-dimensional case 

We have seen that in one dimension, the racemic mix- 
ture naturally tends toward a state where there are 
left and right handed domains. Extending the analy- 
sis is straightforward. Keeping in mind the picture of 
symmetry breaking into different A minima, one can 
envision how the spatial distribution first breaks up 
into different domains. The question is then whether 
these domains are stable or not. This question is eas- 
ily answered by considering a spherical bubble, de- 



scribed by 



d'^A D-l dA 



Q^2 



dr 



(47) 



This has a solution where the RHS vanishes, making 
the solution stable. However, in more than dimen- 
sion {D > 1) this is not the case. Consider the time- 
independent solution in one dimension: a mechan- 
ical analogy is a point mass moving under gravity 
in the inverted potential —V{A) with r playing the 
role of time. The solution that interpolates between 
the two minima corresponds to the point mass rolling 
down from the top of the hill and then stopping at 
the top of the other hill. This is possible because 
in the mechanical analogy, the equations of motion 
include no friction term. In _D > 1 dimensions, the 
situation is different due to the presence of the term 
proportional to {D — l)drA. Analogously to the one- 
dimensional case, this term acts as (time dependent) 
friction. Hence, there is no spherically symmetric 
stable solution in more than one spatial dimension. 
The time dependence will work in the direction as to 
shrink the bubble as is evident from Eq. (|47|) . In a 
more complicated system, where the wall is not a sim- 
ple bubble, this analysis still describes the behavior 
qualitatively: the local curvature determines whether 
a local region tends to shrink or expand. This behav- 
ior is also confirmed numerically. 

Evolution of the enantiomeric excess 

In order to measure the speed of the process, we need 
to consider global quantities that well represent the 
left or right handedness of the whole system. In the 
reduced model we use the enantiomeric excess, as 
defined in Eq. if^ . based on Er = J X{x,t) d^x, 
El = J Y{x,t)d^x. Assuming furthermore S = 1, 
in D dimensions we have r/ = J Ad^x, and so for 
/ = 1 we have 



dr] 
'dt 



J FAix,t) d^x, 



(48) 



because the integral over the laplacian vanishes for 
periodic boundary conditions. Here, L is the size 
of the domain and Z? = 2 or 3 for two and three- 
dimensional systems, respectively. 
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Figure 4: Evolution of ?/ for a two-dimensional sim- 
ulation with k/(L^Ao) = 2 X 10^'', i/Cmicro = IBtt, 
and 1024^ mesh points resolution. The inset shows 
the normalized slope. Note the fours distinct regimes 
with progressively decreasing slope. 



We have calculated the evolution of X{x,t) and 
Y{x,t) through numerical simulations in two and 
three spatial dimension, with and without advection. 
We use the Pencil Code^, which is a modular high- 
order finite-difference code (sixth order in space and 
third order in time) for solving general partial differ- 
ential equations and, in particular, the compressible 
hydrodynamic equations. We usually start with ran- 
dom initial conditions where X and Y are uniformly 
distributed between and 1. The results of the sim- 
ulations are shown in Figsl^HTHl 

In Fig. ^ we show the evolution of 77 for a model 
with large spatial extent; i/^niicro = IStt « 50. Note 
that the evolution proceeds in steps with piecewise 
constant values of 77 (see the inset). The progressively 
decreasing value of f] is associated with a reduction of 
the number M of topologically connected regions in 
the full domain. Looking at Fig. |S1 M decreases from 
A/" = 4 at i/Tdiff = 0.002 to TV = 1 at t/rdis = 0.01, in 
accord with the change of slope of rj{t) see in Fig. ^ 

A blow-up of the little white region in the third 
panel of Fig. [S] [t/T^is = 0.006) is shown in Fig. |H| 
Here, instead of X, we show F4 along with a cross- 
section of the loop. Both inside and outside the ring 
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Figure 5: Plot of X{x,y) at four different times for 
the same simulation as in Fig. 2| White indicates 
X = 1 and dark indicates X = 0. Note that the 
number of disconnected regions decreases from 4 in 
the first plot to 3, 2, and 1. 
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Figure 6: Color/gray scale representation of the 
quantity Fj{{x,y,t) corresponding to the third plot 
in Fig. near the location of the small isolated ring. 
Note the sharp alternations of the sign of F4 just at 
the boundary of the ring. 



Fj[ ~ 0, and F4 ^ only at the location of the 
ring-like front. 

While the fronts shrink and finally disappear, they 
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appear surprisingly robust and the relaxation to the 
final homochiral state is a slow process compared to 
the initial local relaxation to many disconnected re- 
gions with 7] = ±1. So, even though the local re- 
laxation occurs on the time scale Tmicroi the total 
homochiralization process takes place on the much 
longer time scale, Tdiff. 

The piecewise linear scaling of ri{t) is associated 
with the simple front structure shown in Fig. ^ We 
recall that t) is directly proportional to F4. However, 
not only is = outside the front, the positive and 
negative contributions on either side of the front al- 
most cancel. In fact, if the front were straight, the 
two contributions from either side of the front would 
cancel exactly, in which case 77 = 0. The only reason 
the two contributions do not cancel exactly is because 
the inner and outer parts of the front have different 
radii, ri and r2- Denoting the width of each ring by 
Ar, and since their separation is also ~ Ar, we have 
77 (X (r2 — ri)Ar « (Ar)^. Note that (Ar)^ is inde- 
pendent of time. Thus, the value of 77 is independent 
of how complicated the structures are, as long as its 
topology is the same. Furthermore, our simulations 
confirm that (Ar)^ is proportional to k/Xq (see below 
for a quantitative demonstration). 

So far we have always assumed that ki/ks = f = 1. 
If this is not the case, then \/Xq {— 2/ — 1) will be 
different from unity, and so will be \ A\. Instead, A = 
±VX [see Eqs. (43) and (44) of BAHN]. Furthermore, 
the expression for F4 is different from that given in 
Eq. (05, and is now 



Fa=A 



X-A" 
1 + A^' 



(49) 



However, as before, F4 = on either side of the 
front, and a finite contribution comes only from an 
imperfect cancellation from within the front, where 
it is governed by the value of ^ oc a/A- This is also 
confirmed numerically and we find (see Fig. d) 



dry 
dt 



12.57V 



2A^ 
Qfcl 



1/4 



K 
1?' 



(50) 



where we have included the dependence on the num- 
ber of topologically disconnected structures A/", and 
the growth rate Ao has been expressed in terms of Q 




0.0 0.2 0.4 0.6 0.8 
(A/Ao)^/2 



1.0 



Figure 7: Normalized values of 77 as a function of A 
for k/{L^Xo) = 6.3 x 10"^ (diamonds) and 1.3 x 10~^ 
(plus signs). Both groups of data points overlap and 
match the fit formula Eq. H5UI) . In order to reduce 
noise in the date wc have chosen a single ring-like 
structure as initial condition. 
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Figure 8: Normalized values of dij/dt obtained from 
simulations in two dimensions (dots) and three di- 
mensions (circles). In the three-dimensional runs, 
the resolution is varied between 32^^ for k/{L'^Xo) = 
8 X 10-3 and 256^ for k/{L^Xo) = 3 x 10"''. In all 
cases L/^micro = 27r. 



and ks- Note, in particular, that there is no evidence 
for a correction arising possibly from the denomina- 
tor of Eq. (gni). 
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!^/Tturb = 0.1 iJ/Tturb = 0.5 
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Figure 9: Plot of X{x^ y, t) in a two-dimensional sim- 
ulation with turbulence. As in Fig. O kI(1?\q) — 
2 X 10""*, L/^micro = 167r, and 1024^ mesh points 
resolution. Note the narrow elongated structures en- 
hancing significantly the chance of extinguishing one 
or the other life form. 

In three dimensions, the evolution of 7y(i) is no 
longer perfectly linear, but an average value of 77 can 
still be determined. It turns out that these values 
are rather similar to the two-dimensional case; see 
Fig. |S1 In three dimensions the only finite contri- 
bution to Fjs^ comes from a pair of shells with radii 
r\ and r2, so 77 cx (r\ — r^)Ar k. 2r (Ar)^, where 
r is some mean radius that depends on time as the 
shells shrink. This explains the departures from lin- 
ear scaling in the three-dimensional case. In terms of 
K one has i] oc nrjL?. The observed similarity be- 
tween two and three-dimensional cases suggests that 
2r/L = 0(1). 

The main conclusion from these studies is that the 
growth of the enantiomeric excess occurs via diffusion 
and not via front propagation. 

Finally we turn to the case of a finite advection 
velocity u. The evolution equations for X and Y are 
now supplemented by an evolution equation for the 
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Figure 10: Normalized values of dr^/dt obtained 
from simulations with turbulence in two dimensions 
(dots) and three dimensions (circles). In the three- 
dimensional runs, the resolution is varied between 32^ 
for k/(L2Ao) = 3 X 10-3 and 256^ for = 
8 X 10"^. In all cases L/S^^{cra = 27r. 

three components of u and, in this case, for the loga- 
rithmic density. The turbulence is forced in a narrow 
band of wavenumbers that are about 3 times the min- 
imum wavenumber of the box. For the present dis- 
cussion the important quantities are the root mean 
square velocity, Wrms, and the energy carrying scale 
£ {— L/3 in the present case). The details of such 
simulations^ have been discussed elsewhere in more 
detail (e.g., Brandenburg, 2001; Brandenburg et ai, 
2004b). 

In the presence of turbulence, originally detached 
regions with growing right and left-handed polymers 
are driven into close contact, enhancing thereby the 
probability of extinguishing one of the two possible 
file forms; see Fig. O It turns out that 77 still grows 
roughly linearly in time, but now at a dramatically 
enhanced rate; see Fig. ^1 As expected, there is no 
dependence on k and the physically relevant quanti- 
ties are instead Mrms and £. 

It is possible - and indeed customary - to ex- 
press the enhanced turbulent transport in terms 

^Animations of both turbulent and 

purely diffusive solutions can be found at 
jhttp://ww w . nordita . dk/ ~br andenb Tmovi es/ chir al | 
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of a turbulent (eddy) difFusivity. Combining 
77 « (0.3...0.6)urms/^ (from Fig. Uni) and 7) « 
12.5 Ktuih/L'^ [the latter being as a definition of Kturb! 
cf. Eq. (inOJl], we find 

Kturb ~ (0.08.. . 0.15) u^s^, (51) 

where we have used £ — L/i. This is somewhat below 
the canonical estimate of ju^ms^ (Prandtl 1925) for 
turbulent diffusion of a passive scalar, but consistent 
with results of simulations of passive scalar diffusion 
(e.g., Brandenburg et al., 2004b). Thus, the fact that 
the flow has reactive properties does not seem to af- 
fect the turbulent mixing. 

It should also be pointed out that the overall struc- 
ture of X and Y is now rather irregular and no longer 
front-like as in the case of pure diffusion. Thus, the 
fact that we can define a turbulent value of the diffu- 
sivity does not mean that with turbulence the mean 
distributions of X and Y can accurately be described 
by a diffusion term. 

Discussion 

Our studies have shown that, without advection, the 
enantiomeric excess grows toward homochirality with 
the speed t) ^ k/L^, where L is the physical size of 
the system; see Eq. (|50|) . This can be much slower 
than the propagation of epidemics which proceeds on 
the faster time scale Tfront — -^/"f^ftont- For the Earth, 
this would yield a time scale for homochiralization 
of T 10^"^ s, i.e. comparable to the age of the uni- 
verse. So the diffusive process is too slow in creating 
large homogeneous areas of one chirality, even though 
within a local region, a single handed state is created 
on the time scale A^^^, which could be on the order of 
days [see the discussion after Eq. (|24|l ] , or even years 
(as in the context of the spreading of epidemics). 

Turbulence, however, changes the situation drasti- 
cally: the speed toward homochirality is effectively 
77 ~ ttrms/ij which, assuming a root mean square 
flow speed of say 1 cm/s, leads to time scales of a 
hundred years, in rough agreement with results from 
ocean mixing studies where horizontal eddy diffusiv- 
ity is found to be of the order of lO'^m^/s (Stewart 
2003, Chap. 8.5). 



Comparisons to present day oceanographic stud- 
ies must be done with care, however, as there are 
a number of factors that can affect both horizontal 
and vertical mixing such as the higher salinity in the 
early oceans (Knauth 1998) and most importantly 
the possible lack of continents on the early Earth. 
The role of continents is manifold as they aid in the 
vertical mixing and are possibly needed in creating 
large scale flows (Garrett, 2003). Hence, the condi- 
tions and mixing in the early oceans may have been 
quite different from the complex system of currents 
what we observe today. In addition, frequent impacts 
on the early Earth may not only have led to increased 
mixing, but also to global extinction of life more than 
once (Wilde et al, 2001). 

Although it is impossible to know with any degree 
of certainty the amount of mixing in the early ocean, 
the time scales may have been rather long. If this 
is the case, the question of homochirality would re- 
main open. Thus, it might therefore be interesting 
to investigate in more detail whether the chirality of 
biomarkers from early marine sediments can be de- 
termined. An obvious difficulty with this proposal is 
the fact that dead matter tends to racemize sponta- 
neously on time scales that are typically shorter than 
a million years (e.g., Bada, 1995). 

Acknowledgments 

The Danish Center for Scientific Computing is ac- 
knowledged for granting time on the Linux cluster in 
Odense (Horseshoe). 

References 

Bada, J. L. (1995) Origins of homochirality. Nature 
374, 594-595. 

Brandenburg, A. (2001) The inverse cascade and non- 
linear alpha-effect in simulations of isotropic 
helical hydromagnetic turbulence, Astrophys. 
J. 550, 824-840. 

Brandenburg, A., Andersen, A., Nilsson, M., 
and Hofner, S. (2004a) Homochiral 
growth through enantiomeric cross- 



13 



inhibition, Orig. Life Evol. Biosph. (in 
press). Preprints available online at: 
jhttp : //arX iv . org/abs/q-bio/0401036 
(BAHN). 

Brandenburg, A., Kapyla, P., and Mohammed, 
A. (2004b) Non-Fickian diffusion and tau- 
approximation from numerical turbulence, 
Phys. Fluids 16, 1020-1027. 

Cleaves, H. J. and Chalmers, J. H. (2004) Ex- 
tremophiles may be irrelevant to the origin of 
Hfe, Astrobiol. 4, 1-9. 

Collins, L. J., Poole, A. M., and Penny, D. (2003) 
Using ancestral sequences to uncover potential 
gene homologues, Appl. Bioinform. 2, S85- 
S95. 

Cotterill, R. (2002) Biophysics: an introduction, J. 
Wiley & Sons, Chichester. 

Fisher, R. A. (1937) The wave of advance of advan- 
tageous genes, Ann. Eugenics 7, 353-369. 

Frank, F. (1953) On Spontaneous Asymmetric Syn- 
thesis, Biochim. Biophys. Acta 11, 459-464. 

Garrett, C. (2003) Internal tides and ocean mixing. 
Science 301, 1858-1859. 

Haken, H. (1983) Synergetics - An Introduction, 
Springer: Berlin. 

Joyce, G. F., Visser, G. M., van Boeckel, C. A. A., 
van Boom, J. H., Orgel, L. E., and Westrenen, 
J. (1984) Chiral selection in poly(C)-directed 
synthesis of oligo(G), Nature 310, 602-603. 

Kaliappan, P. (1984) An exact solution for travelling 
waves of ut = Duxx + u — u'^, Physica D 11, 
368-374. 

Knauth, L. P. (1998) Salinity history of the Earth's 
early ocean. Nature 395, 554-555. 

Kozlov, I. A., Pitsch, S., and Orgel, L. E. 
(1998) Oligomerization of activated D- and L- 
guanosine mononucleotides on templates con- 
taining D- and L-deoxycy tidy late residues, 
Proc. Natl. Acad. Sci. 95, 13448-13452. 

Martin, W. and Russell, M. J. (2003) On the origins 
of cells: a hypothesis for the evolutionary tran- 
sitions from abiotic geochemistry to chemoau- 
totrophic prokaryotes, and from prokaryotes 
to nucleated cells, Phil. Trans. Roy. Soc. B 



358, 59-83. 

Martinez, J. S., Zhang, G. P., Hoh, P. D., Jung, 
H.-T., Carrano, C. J., Haygood, M. G., 
Butler, A. (2000) Self-assembling amphiphilic 
siderophores from marine bacteria. Science 
287, 1245-1247. 

Mendez, V., Fort, J., Rotstein, H. G., and Fedotov, 
S. (2003) Speed of reaction-diffusion fronts in 
spatially heterogeneous media, Phys. Rev. E 
68, 041105. 

Milner- White, J. E. and RusseU, M. J. (2004) Sites 
for phosphates and iron-sulfur thiolates in the 
first membranes: 3 to 6 residue anion-binding 
motifs (nests), Orig. Life Evol. Biosph. (in 
press). 

Mojzsis, S. J., Arrhenius, G., McKeegan, K. D., Har- 
rison, T. M., Nutman, A. P., Friend, C. R. L. 
(1996) Evidence for fife on Earth before 3,800 
million years ago. Nature 384, 55-59. 

Murray, J. D. (1993) Mathematical Biology, Springer, 
Berlin. 

Nelson, K. E., Levy, M., and Miller, S. L. (2000) 
Peptide Nucleic Acids rather than RNA may 
have been the first genetic molecule, Proc. 
Nat. Acad. Sci. U.S.A. 97, 3868-3871. 

Nielsen, P. E. (1993) Peptide nucleic acid (PNA): A 
model structure for the primordial genetic ma- 
terial, Orig. Life Evol. Biosph. 23, 323-327. 

Pal, D., Siihnel, J., and Weiss, M. S. (2002) New 
principles of protein structures: nests, eggs- 
and what next? Angew. Chem. Int. Ed. 41, 
4663-4665. 

Pooga, M., Land, T., Bartfai, T., Langel, U. (2001) 
PNA oligomers as tools for specific modulation 
of gene expression, Biomol. Eng. 17, 183-192. 

Prandtl, L. (1925) Bericht iiber Untersuchungen 
zur ausgebildeten Turbulenz, Zeitschr. Ange- 
wandt. Math. Mech. 5, 136-139. 

Rasmussen, S., Chen, L., Nilsson, M., and Abe, S. 
(2003) Bridging nonliving and living matter, 
Artif. Life 9, 269-316. 

Russell, M. J. and HaU, A. J. (1997) The emergence 
of life from iron monosulphide bubbles at a 
submarine hydrothermal redox and pH front, 
J. Geol. Soc. Lond. 154, 377-402. 



14 



Sandars, P. G. H. (2003) A toy model for the gener- 
ation of homochirality during polymerization, 
Orig. Life Evol. Biosph. 33, 575-587. 

Saito, Y. and Hyuga, H. (2004a) Complete homochi- 
rality induced by the nonlinear autocatalysis 
and recycling, J. Phys. Soc. Jap. 73, 33-35. 
(SH) 

Saito, Y. and Hyuga, H. (2004b) Homochirality 
proliferation in space, J. Phys. Soc. Jap. 
(submitted). Preprints available online at: 
'http : //arXiv . org/abs/physics/040506^ 

Schopf, J. W. (1993) Microfossils of the early archean 
apex chert - new evidence of the antiquity of 
hfe. Science 260, 640-646. 

Stewart, R. H. (2003) Introduc- 

tion to Physical Oceanography, 

'http : / / oceanworld . tamu . edu/ resources/oc ng_textbook/ contents . html| 

Tedeschi, T., Corradini, R., Marchelli, R., Pushl, A., 
and Nielsen, P. E. (2002) Racemization of chi- 
ral PNAs during solid-phase synthesis: effect 
of the coupling conditions on enantiomeric pu- 
rity. Tetrahedron: Asymmetry 13, 1629-1636. 

Ward, P. D. and Brownlee, D. (2000) Rare Earth: 
why complex life is uncommon in the Universe, 
Copernicus, Springer: New York. 

Watson, J. D. and Milner- White, E. J. (2002) A novel 
main-chain anion-binding site in proteins: the 
nest. A particular combination of ^,"0 val- 
ues in successive residues gives rise to anion- 
binding sites that occur commonly and are 
found often at functionally important regions, 
J. Mol. Biol. 315, 187-198. 

Wattis, J. A. D. and Coveney, P. V. (2004) 
Symmetry-breaking in chiral polymer- 
ization, Orig. Life Evol. Biosph. (in 
press). Preprints available online at: 
|http : //ar Xiv . org/abs/physics/040209i| 

Wilde, S. A., Valley, J. W., Peck, W. H., and Gra- 
ham, C. M. (2001) Evidence from detrital zir- 
cons for the existence of continental crust and 
oceans on the Earth 4.4 Gyr ago. Nature 409, 
175-178. 

Id: paper. tex,v 1.90 2004/07/05 20:08:30 hrandenb 
Exp 



15 



